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Abstract 

These notes describe in some detail the continuum limit behavior of a very simple car following 
traffic flow model. The formation and behavior of shock waves is described. This model is the 
one solved by a set of MatLab scripts in the Athena 18311-Toolkit at MIT, which illustrate the 
phenomena described here. These are the scripts whose names end with the acronym CFSM. 



Contents 

1 The Model. Nondimensionalization. 2 

2 Continuum Limit of Model. 4 

3 Numerical Issues. Stiffness of the equations. 9 

4 Examples. 13 

5 Notes on the MatLab script quadCFSM. 15 

List of Figures 

1.1 Typical flow and velocity functions 2 

5.1 Typical initial conditions for script quadCFSM 15 

5.2 Solution after shocks form 17 

5.3 Approximation of the initial wave velocity data 19 

5.4 Comparison of exact and approximate solutions 21 

*MIT, Department of Mathematics, room 2-337, Cambridge, MA 02139. 

1 



Simple Traffic Flow Model. MIT, Friday March 26, 1999 — Rosales. 

1 The Model. Nondimensionalization. 

Consider a line of cars on a road, with car n located at x n = x n (t), moving at speed u r 



~dJ 



Measure distance x along the road in the same direction the cars move (so the car velocities u n are all 
non-negative). Number the cars so that {x n } is an increasing sequence (x n +i — %n) > car length > 
(identify x n with the location of (say) the front end of the car). 

Remark 1.1 We use tildes over the variable symbols to indicate that we are dealing with dimen- 
sional variables. We will use the same symbols (without tildes) to denote the nondimensional ver- 
sions of the same variables, when we nondimensionalize the equations later on. 

Assume now that the drivers follow some rule (such as the ones proposed in problems 61.1 and 61.2 
of Haberman's book) prescribing the car speed as a function of the distance to the next car 1 . That 
is, assume that there is some function (the car velocity function) 



U = U(p) 



such that 



U r 



U (p n ) , where p m 



1 



(1.1) 



Here we identify p n with the car density at the position of the n th car. We also introduce the 



Typical function u = U(p) 




Typical function q = Q(p) 




Figure 1.1: Typical forms for the car velocity and flow functions. 



1 As long as the situation is not changing too rapidly, this is not unreasonable. Note that in this model we will, 
implicitly, deal with all the cars as if they were equal copies of each other — all the cars obey exactly the same rules. 
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notation h n = (l/p„) = x n +\ — x n for the car separation. Typical shapes for the car velocity U and 
the car flow Q = pU (both functions of p) are shown in figure 1.1. 
Then the model is given by the following set of coupled ODE's 

1r = ^' (L2) 

where the velocities and positions are related by equation (1.1). To complete the model we need to 
give a boundary condition. For example, if there are N cars, then the velocity un of the car at 
the head of the group would have to be prescribed (since (1.1) cannot be used for n = N). 

Typical values 2 for the various constants involved are as follows: jamming density pj = 160 cpm | , 



road capacity | q m = 1600 cph | , density at road capacity | p m = 80 cpm | and car velocity at road 
20 mp/Tj . If we now assume a length scale L — characterizing a typical length over 



capacity 



which the traffic density changes significantly — we can nondimensionalize as follows: 

x n = Lx n , p~ n = pjp n , u n = —u n , U(p) = —U(p) , Q(p) = q m Q(p) and t = — t . 

PJ PJ Qm 

(1.3) 

Then the equations take the form 

u n = U(p n ) and p n = , (1.4) 



at x n +\ x n 

where e = \j{Lpj) is a small nondimensional number — with the values above and with L a large 
fraction of a mile, we get 



e = O(10~ 2 ) 



Note also that the nondimensional versions of the car 
velocity and car flow functions have the same forms as the dimensional ones; but with the jamming 
density and road capacity set to one. 

Remark 1.2 In the nondimensionalization above, the choice of L was left a bit ambiguous. While 
the other parameters follow from actual measurements and are pretty fixed (for a given road), the 
length scale is more flexible and depends on the particular solution of the equations one is looking at. 
On the other hand, while L can (at least in principle) be arbitrarily large 3 , there is an approximate 
minimum size | L min | it can have. Consider the way we defined L, which requires (in particular) 
that we be able to distinguish a length scale. Thus, think of a typical perturbation to the density 



2 See sections 62 and 63 in Haberman's book 

3 Consider the example where all the cars are equally spaced, so that L = oo. 
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p along the road — say, a hump or a sinusoidal up and down. This perturbation is "marked" by 
a discrete set of points (the car positions) and needs a minimum number of them before it can be 
clearly identified — a reasonable number 4 being about twenty or so. With the kinds of car densities 
implied by p m , we see that | L min | cannot be much shorter than about a quarter of a mile 
(best we can have, at near jamming density, is about an eighth of a mile). Thus our assumption 
above (where we took L a large fraction of a mile) is quite reasonable. This, in addition to 



0(10 2 ) , yields (in the nondimensionalization above in (1.3)) a time scale 



order of a few minutes — 6 min for L a full mile and 3 min for half a mile. 



Lpj 



in the 



Remark 1.3 Continuing with the issue of the length scale L: we also need (of course) that a length 
scale exist! The cars could be randomly distributed on the road, in which case there would not be 
much of a length scale to be identified. However, this a situation that cannot persist: think of the 
example of three cars, with the first two close and the third far behind. Then the third car would 
end up moving faster than the second and the two distances would tend to even out. If, on the other 
hand, it is the second and third car that are closer, then the second car would move faster than the 
third and (again) the distances would tend to even out. In general, there is a tendency for the cars to 
settle down to situations where the car separations do not vary very rapidly, except for a few isolated 
places where "jumps" occur. This process is illustrated by the MatLab script randCFSM in the 
Athena 18311-Toolkit, which solves the equations in this model with random initial separations 
between the cars. We will come back to these issues in remarks 3.1 and 3.3. 



2 Continuum Limit of Model. 

In remark 1.2 we considered the question of what is the minimum size the length scale L — used in 
equation (1.3) to nondimensionalize lengths — can have. In this section we will consider a somewhat 
opposite situation, where L becomes larger and larger. Equivalently: 

consider the limit: |e — [ 0| in the model equations (1.4). (2.1) 

Of course, this is only an "ideal" limit we are taking. In practice e is fixed. However, since e is 
small, we expect the limit will give us useful information regarding the behavior of the model (1.4). 

4 Think of how many points per wavelength are needed to have a reasonable drawing of a sine wave. 
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Remark 2.1 Let N L be the typical number of cars that make the features (bumps, whatever) in the 
traffic density that were used (earlier, in equation (1.3)) to determine the length scale L. These two 



quantities are related by an equation of the form N L = p x L , where p x is some average density 5 



(which cannot be too different from, say, p m ). Thus, letting e — > amounts to both making L and 

N L large, since e = - — = ^* . In fact, note that 
Lpj pjN L 



e = 0(N£ 1 ) - since ^ = 0(1). 

Pj 



Because of the way the equations were nondimensionalized, we see that: 

• The separation between cars satisfies 

h n = x n+1 -x n = 0(e) . (2.2) 

This follows because in equation (1.4) p n = 0(1). 

• Significant variations in car separation (i.e., in p n ) occur over 0(1) distances. Thus, it is 
reasonable to assume that there is a function p = p(x, t) such that 

pn = p(x n ,t). (2.3) 

We expect p to be reasonably nice and (generally) have 0(1) partial derivatives — and — . 

We now rewrite the equations for the model (1.4) in terms of the densities rather than the car 
positions. Thus we have 



^r.pn = -e 1 pn ( u n+i ~ u n ) or, equivalently: ^-p n = -p n 1^+1 ^ J ? (2.4) 

Cut CVt \ Xn-\-i Xyi J 



where the densities, velocities and positions are related, in the usual way, by u n = U(p n ) and 
Pn = e/(^n+i — x n ). Again, in addition to initial conditions a boundary condition is needed. For 
example, if there are N cars, then velocity (or the density) u N of the leading car would be required. 

From equations (2.2 - 2.3) it is clear that the expression [ n+l J can be replaced by —(x n , t) 

in the limit given by (2.1), where u = u(x,t) = U(p(x,t)). Furthermore, from (2.3) it follows that 

—p n = ( — I- I (x n ,t) — using the chain rule. Substituting all this into equation (2.4) we see 
at \ot ox ) 

that p above in (2.3) must satisfy the PDE 

= ^ + ^ = ^ + c^ (2.5) 
at ox at ox 



'For example: we used p* = p m in remark 1.2 (with Nl « 20) to determine L ri 
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in the limit e — > 0, where q = pu = pU(p) = Q(p) and c = c(p) = -j-. Thus we obtain the same 
continuous traffic flow model that was developed in the lectures (see the lecture notes or the book 
by Haberman) using a phenomenological approach and conservation of cars. 

An interesting point arises now. The solution of the PDE (2.5) (by characteristics, say) generally 
breaks down after a finite time. That is, infinite derivatives and multiple values develop after some 
critical time — even if the initial data are smooth. On the other hand, it is quite clear that the 
model (1.4) — equivalently (2.4) — cannot develop anything even resembling multiple values. In 
fact, there is no breakdown either: provided the initial values for the car positions are such that 
the densities all satisfy < p n (0) < 1, then the solution will exist for all times and the bound 
< p n (t) < 1 will be satisfied. The proof of this is rather easy: (i) The density can go to zero only 
if the distance between cars goes to infinity, but this cannot happen because the car velocities are 
bounded, (ii) Neither can the density go beyond one, for as soon as p n reaches one, the n th car will 
stop, while the (n + l) th car will be moving at a non-negative velocity, (iii) Thus, the condition 
< p n < 1 will be preserved, (iv) This is enough to guarantee a solution for all times, for a solution 
can cease to exist only if it either "blows up" or if it reaches a singularity in the equations. However, 
as long as < p n < 1, neither of these two things can happen. 

Note 2.1 Notice that the argument in (ii) above shows that a density of one can be maintained 
only if the density is identically one from some car on forward. Else a decrease in density will 
propagate backward through the cars, as the cars where p n = 1 will not move. On the other hand, if 
the density is one from some point on, then a "wave" carrying a value of one will move backward 
through the road, as cars move into the ones that are stopped ahead in the line. 

The question is now: what happens in the limit (2.1) beyond the time where the solution 
of the PDE (2.5) breaks down? This question is addressed by the MatLab script quadCFSM 
in the Athena 18311-Toolkit — see section 5 here. This script solves the model equations (1.4) 
— with initial conditions that correspond to a smooth positive hump for (2.5) — in the limit (2.1). 
Actually not "in the limit", but for e small enough that one can see what will happen when e — > 0. 
With the initial conditions used by quadCFSM the solution to (2.5) breaks down in a finite time. 
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On the other hand, what the numerical experiments show is: 

1. As long as the solution of (2.5) behaves nicely, it does approximate quite well 
the behavior of the solution of (1.4). 

2. The solution of (2.5) exhibits breakdown with formation of infinities in the 
derivatives in the regions where the density p is increasing with x. In these 
regions, the solution of (1.4) also shows progressive steepening of the density 
profile. However, rather than "topple over" and develop multiple values 
(as happens with the solution by characteristics of (2.5)), the solution of f (2.6) 
(1.4) develops a very sharp transition — just a few cars wide — from one 
value of p to a bigger one. In effect, the function p in (2.3) develops a 
discontinuity that stops multiple values from arising. 

3. Other than the phenomena described in the prior item, the e — > limit of 
(1.4) is still described by (2.5) — that is, away from the discontinuities in 
the density, (2.5) applies. 

Thus the proper description of the limit in (2.1) is still (2.5), but we must add discontinuities 
(across which the density increases) in the solution to avoid the formation of multiple values. 
These discontinuities are called SHOCKS and cannot be placed arbitrarily, since: 

• Shocks must move so that cars are conserved. If x = xg{t) is the shock position, then 



d 

dt Xs 



l_0\ 



(2.7) 



where we use the notation [ ] to denote the jump in a function across a discontinuity. This 
condition is called the Rankine— Hugoniot jump condition. 
The (so called) entropy condition must hold 



across shocks the density increases. 



(2.8) 



In terms of the characteristic curves for equation (2.5), this means that the curves converge 
into the shock — and terminate there. Thus the shock path acts as a "cut" in space- 
time, where the characteristic curves end. This prevents their crossing and the formation 
of multiply valued regions in the solution. 
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It can be shown that these two conditions are enough to uniquely determine the solution of (2.5), 
now for all times and without any multiple values arising. Thus, this "augmented" model (i.e. 
equations (2.5), plus discontinuities governed by (2.7) and (2.8)) is the result of (2.1). 

Remark 2.2 The condition (2.8) is very important. No discontinuous transitions are devel- 
oped by (1-4) that are not associated with an increase in the density. This is very clear intuitively; 
when the density is decreasing the cars move faster the further ahead they are in the line and no 
steepening tendency arises (exactly the opposite occurs). It is only when the density increases that 
sharp transitions are generated and maintained. 

We now examine the derivation leading to equation (2.5) and ask: what did we miss that would 
explain the behavior in (2.6)? The answer has to do with the assumption right below (2.3) that 
p has 0(1) partial derivatives — which it obviously does not. Thus there will be extra contributions 
(that we neglected) near shocks to equation (2.5). Specifically, consider the step where we replaced 

(— J by —(x n , t). In a more precise calculation (to estimate the error made by the substi- 

tution) we expand u n+ \ in a Taylor series centered at x n . That is u n+ \ = u n + v$h n + -v$h 2 n + . . ., 

At 

where h n — x n +\ — x n and we use the notation u$ = 7— (x n ,t). Thus 

0X-* 

fer^) = I'*--" + ^ + s« + = <> + ^ + 4 e2 "" 1 + • 

where we used (1.4) to replace h n by p n . In this expression we can neglect the higher order terms 
(as we did when obtaining (2.5)) only if we can argue that e 7-1 ^^ is small for j > 1. But, when 
the partial derivatives of p — thus those of it — are not bounded, we cannot do this. Then all these 
extra terms will become important (an cannot be neglected, as we did in (2.5)) near shocks. 

To have an idea of what the effect of these extra terms is on the behavior of the solution, it is 
enough to just keep one extra term and see how this changes (2.5)). Using u = U(p) to replace 
derivatives of u by derivatives of p, this yields the equation 



where 



dU 
dp 



(notice that v is a POSITIVE function of p). Thus a (small) amount of 



diffusion is added to equation (2.5). As long as the derivatives are bounded, the effects of this 
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diffusion can be neglected. But, when the density profile steepens, it becomes important and begins 
to "fight" the steepening — this is what diffusion does. Eventually a balance between the two effects 
(diffusion and the tendency to steepen) is achieved, within a narrow region of high derivatives. 

It is easy to estimate the width the balanced region in the prior paragraph should have, as follows 
(this will be the shock width). Let this width be ws- Then, while p will remain 0(1) near the 
shock, each derivative will be larger than the prior one by a factor Wg 1 . Thus, in equation (2.9), 
the left hand side will have size Wg 1 while the right hand side has size ewg 2 . Clearly, for balance we 
need wg = e- Since e is also the order of magnitude of the car separation, this predicts a shock 
width of a few cars, which agrees well with the numerical results in (2.6). 



3 Numerical Issues. Stiffness of the equations. 

We now go back to the discrete equations and perform an analysis to see what sort of time scales 
are involved in their behavior. This is important for many reasons, some of which we will explain 
later on. In particular: in any numerical calculation we must make sure that all times scales are 
handled properly, even if they are not immediately apparent in the solution — the precise meaning 
of this last rather strange statement will be clarified below in remark 3.2 (second point). 
Consider a situation where the car densities deviate slightly from some constant state p*. Thus: 

Pn = P* + Sn , (3-1) 

where the perturbations S n to the density are assumed small. Substituting this formula in the 
equations for the model (use the first form in equation (2.4), which involves only the densities p n 
and the velocities u n = U(p n )) and neglecting higher order terms in the perturbations, we obtain: 

^ = t~ x v*p\ (S n+1 - 6 n ) , (3.2) 

where v x = v(p*) and v is as in (2.9). This last equation is linear and can be solved using eigenmodes. 
Specifically, the general solution is a linear combination 6 of the fundamental modes: 

S n = e ( ikn + Gt ) , with a = e- l v*pl (e ik - l) = e"Vp* {(cos(fe) - 1) + isin(k)} , (3.3) 
6 Notice that this is the same type of solution used in the von Neumann stability analysis of numerical schemes. 
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where — n < k < n (these solutions are periodic in the wavenumber k, since the exponential is 

sampled only at integer values). The values of a determine the time scales involved in the solution. 

We note that all the er's have negative real parts, so that all these solutions decay (i.e. the constant 

state p x is a stable solution of (2.4)). In fact, the shorter the wavelength A = — , the faster the 

k 

decay rate. The maximum decay rate corresponds to solutions that oscillate with a wavelength of 
two car separations (k = ±7r), with a = —2e~ l v^p\. This corresponds to a time scale 

Remark 3.1 We note that r m is a very short time. As pointed out in remark 1.2, 0(1) times 
in the nondimensional equations typically correspond to a few minutes in dimensional units. Since 
e = O(10~ 2 ) ; we see that r m corresponds to a dimensional time scale that must be measured 
in seconds! Now we ask and answer the questions: What exactly is the meaning of the time 
scale r TO ? What role does it play in the time evolution of the equations? For this we go 
back to the point made in remark 1.3. It is quite clear that r m is precisely the time scale over 
which rapid variations in the car separations are "wiped out" by the time evolution of 
the model. This is the process illustrated by the MatLab script randCFSM. After these variations 
are eliminated, this time scale plays no role, except to the extent that it keeps eliminating any such 
small variations that might arise due to "external" perturbations. 

Remark 3.2 The last statement in the prior remark appears innocuous, but it is actually not. 
What do we mean here by "external" perturbations? 

• First: the equations (1-4) are a pretty crude model for traffic flow; it is pretty unrealistic to 
assume that the drivers respond only to the distance to the car right ahead (and then that they 
can adjust their car velocity instantaneously to the prescribed u). We are using this model 
only as a simple example to illustrate some of the phenomena involved. However, even if we 
were to set up an ideal situation, it would still be an approximation. Thus, all the neglected 
"little" things that the model ignores would constantly introduce changes (perturbations) into 
the solution. In addition one would still have to consider truly external perturbations, such as 
a new car added (or one gone) to the line. Note that it is important that a mathematical model 
be "stable" to such perturbations, else it is worthless (as the neglected effects would be able to 
completely change the nature of the solution). On this last account (at least), the model (1-4) 
behaves the right way. 
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• Second: another (very important) source of "external" perturbations arises when solving 
the equations numerically. This is because any numerical scheme will, necessarily, in- 
volve approximations — which will introduce errors into the solution. These errors better 
not grow, else disaster will strike. Now, the exact equation here would very quickly dissipate 
them, but this need not be so with a numerical scheme if one is not careful. Precisely because 
the equations being approximated are so forceful about dissipating errors, naive numerical 
approximations can easily over do the effect and end up amplifying the pertur- 
bations! A simple example of this is provided by the equation | y = — ay~\ with a large 
and positive. The solutions of this simple equation decay very fast to zero. But, approxi- 
mate the equation by the naive forward Euler scheme: 



Vn+1 = V n - ay n At, 



where 



y n = y(nAt). Then y n = (1 — aAt) n y and, unless At < 2/ a, the numerical solution 

blows up! Thus, to get this scheme to behave properly one needs to take a time step which is 
as short as the time scale of decay. For the equations given by (1-4) this would mean a time 
step as short as r m , which is disastrous! That is, we would be forced to resolve time scales 
in the order of seconds (or fractions), while in fact the phenomena we are really interested in 
following take place over minutes or even hours. In fact, nothing happens over seconds, we 
have to keep such a small time step just so the numerical scheme does not go unstable. 

Problems that present short time scales that are irrelevant to the solutions one is trying to 
compute (but arise because small deviations from these solutions are very quickly "squashed" 
by the governing equations) are called numerically STIFF and require special care. Naive 
approximations invariably lead to very inefficient codes, requiring unrealistically small time 
steps. We will not go into these problems here, but you should be aware of their existence. 

Remark 3.3 Actually, as pointed out earlier, all the a's in (3.3) have negative real parts, so all 
the scales decay. Thus, if we wait long enough, not just the short wavelength (a few car distances 
long) variations will vanish, but the long ones as well. Although this conclusion is based on the 
linearized analysis in (3.1 -3.3) and thus is valid only for small perturbations from a constant, it 
is actually true for the whole set of equations (as a bit of numerical experimentation will quickly 
show). This then makes it necessary that we revisit the statements made in remark 1.3, and state 
them in a more precise way. The "natural" state for the model is to go into a situation where the 
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12 



length scale satisfies L = oo. Any length scale present in the initial data will eventually be wiped 
out (unless it keeps on being re-introduced by external perturbations). But the large scales have 



while short scale variations will be quickly dampened (and will become irrelevant), longer scales will 
remain for "reasonable" times. Thus, we are back to being rather vague about the meaning of the 
space scale L. Basically, we have to argue phenomenologically: it is produced by processes that are 
very complicated and are not included in the model. At the level of simplicity of this model, there 
is not much more that we can do about it. We must take this scale L as an external input, on the 
same footing as other quantities such as pj, etc. The value for L min computed earlier gives an idea 
of what is reasonable (i.e. anything larger than L min and smaller than the length of the road), but 
this is about as much as we will be able to say here. 

Remark 3.4 One can make an interesting observation regarding the size of r m . It is clear that 
the model (1-4) does not allow accidents (car collisions). These would require (at the very least) 
that p n >l somewhere, sometime. But we showed earlier (see the paragraph above the note 2.1) 
that the equations will not let this happen. The time r m is closely associated with the mechanism 
that prevents this from happening. Now notice that real accidents happen 7 (even when the drivers 
attempt to follow the recommended rules of separation between cars versus speed) because of effects 
we have not considered in the model, such as: (a) human reaction time, (b) cars cannot accelerate 
or stop instantaneously, etc. Unlike the mechanism behind the time r m , which is stabilizing, these 
other effects destabilize. The interesting fact is that the time scales associated with them are about 
the same as those given by r m . But perhaps this is not too surprising, if one postulates a tendency 
to "push the envelope" in terms of safety. That is: drivers will drive as fast and as close to the 
next car as it is "reasonably" safe, where this means that the stabilizing effects will be kept at a 
"multiple" of the de-stabilizing ones — but not too large a multiple! 

Finally, just for completeness, we end this section by showing how the linear perturbation analysis 
in (3.1 - 3.3) looks like in terms of the car positions. In this case we have 



decay times much longer than r, 



in 



since the real part of a behaves like k 2 for k small. Thus, 



e 

x = x + n h uj + y,, 

P* 



(3.4) 



Let us exclude here such things as the drivers falling asleep, etc. 



Simple Traffic Flow Model. MIT, Friday March 26, 1999 — Rosales. 

where u* = U(p*) and y n is small. Then 

e e e e 

Vn+l Vn H ~^^n t 

P* Pn P* Pi. 

where we used (3.1) and neglected quadratic terms in the perturbations. Thus 

pl 

$n = J (Vn+l ~ Vn) ■ 

Since u n = U(p x + 5 n ) = u x — v x 5 n , we then have (using (3.4) 

dy n _ 1 

~TT — o [Vn+l — Vn) ? 
dt I T m 

which (of course) is the same equation the S^s satisfy! 

4 Examples. 

In this section we consider examples of choices for the velocity U(p) and flow Q(p) functions. We 
stress that these are just qualitative examples, not actual fits to measured data (which need not 
give simple formulas). Thus one must be careful about not drawing too many conclusions from 
them, specially of the too precise quantitative type. 

Example 4.1 The simplest example is that of a quadratic flow function, 

Q = p (pj - p) and U = (pj - p) . 
Pj Pj 

This yields p m = - pj, u m = — — and a maximum car velocity u max = 2u m . These numbers are 
2 pj 

compatible with the typical values given earlier above equation (1.3), except that the maximum car 
velocity seems a bit low (though not out of range). Then again, the typical values given are from 
measurements in the NYC Lincoln tunnel in the 1950's (where, perhaps, a maximum car speed 
of 40 mph was reasonable). In general these numbers are meant only as ballpark figures. After 
nondimensionalization, we have the forms 

Q = 4p(l-p), U = A{l-p) and c = 4(1 - 2p) . (4.1) 

In this case the shock speed in (2.7) is the average of the characteristic speed c across the shock and 
v = 4 in (2.9). In particular r m = —r. 
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(3.5) 
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Example 4.2 Another simple example follows from the rule stating: for each unit v r of some speed 

(v r = 10 mph is typical) the separation to the next car should increase by at least one car length I. If we 

u n 1 
apply this rule exactly, then — I + 1 = x n+ \ — x n = — . From this and the speed limit, we obtain 

V r " Pn 

u = U(p) = min (u max , v r PJ „ and Q = pU = min (pu max , v r (pj - p)) , 

where pj = l~ x . This yields 

_ _ Vr PJ , _ _ Vr PJ U max 

v^m Vj max , pm ana, Q m p m u m 

Vr i Vi max V r -]- U max 

With u max = 50 mph, v r = 10 mph and pj = 160 cpm this yields p m ~ 27 cpm and q m ~ 1330 cph 
— not altogether unreasonable numbers. One point though is that pj = 160 cpm corresponds to 
I = 33 ft, which is a tad too long. The reason for this is that the cars stop when the distance to 
the next car is bigger than zero (not zero, as this rule would have). Thus, if one uses pj = l~ x with 
an actual car length, too high a jamming density results — so we use a car length that is about 
twice actual to compute pj. In other words, this rule is rather unrealistic for low velocities. The 
implementation of the speed limit is also rather crude and gives the strange feature of a corner (at 
the maximum) in the flow profile Q = Q(p). After nondimensionalizing, we have 

Q = min (£ , i^£) , U= m,n ( - , -IziU and c = L^+iMpj , (4 .2) 
\a 1 — a) \a p (1 — a) ) I a (1 — a) 

where < a = — < 1. Note the strange feature of a piece-wise constant wave speed 

Pj 

c. Thus, in the continuum limit, the parts of the density profile with p > a move (backwards) at 
constant? speed (a — Similarly, the parts with p < a move (forward) with speed aT 1 . Shocks 

will arise where these two kinds of behaviors "collide" and will move at the speed given by (2.7), 
with a jump in p (as x increases) from p < a to p > a. This is pretty strange behavior! This case 
is implemented by strangeCFSM in the Athena MatLab 18311-Toolkit. 

Finally, note that an alternative formulation of the rule in this example is that the time it would 

take a car to cover the distance to the next car should be a given fixed At. It is easy to see that 

£ 1 

the correspondence is v r = — — since this rule simply states that uAt = £. In particular. 

At p 

v r = 10 mph and I = 16 ft correspond to At = 1.1 sec, since a mile is 5280 ft. 



^Therefore: no wave shape deformation. 
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5 Notes on the MatLab script quadCFSM. 

The MatLab script quadCFSM in the Athena 18311-Toolkit solves the equations in (1.4) us- 
ing the quadratic flow function (4.1) in example 4.1. A finite number of cars N is used, with 
X\ < x 2 < ■ ■ ■ < xn and the density < pN < 1 at the leading car given and constant 9 . The ini- 
tial conditions are such that (see figure 5.1) x N (0) = 0, Xi(0) < —it and 



Pn(0) = Pn + (1 - pN)r{x n ) for 1 < n < N 



(5.1) 



with r = r(x) a symmetric positive "hump" in — tt < x < 0, r(— 7r) = 1 and r = outside [— 7r,0]. 

At 



Typical initial conditions: p n = p(x n ) 




Figure 5.1: Typical initial conditions quadCFSM. 



Remark 5.1 The cars are placed so that x p (0) = —n for some p > 1 (thus N h = N + 1 — p is the 
number of cars in the hump, with 1 < Nh < N). Then, from (1-4), we must have 

N-l e 

7T = X N (0) - X p (0) = . 

n=p Pn 



3 The leading car velocity is then also constant un = 4(1 — pjy)- 
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This equation determines the value of e in terms of the number of cars in the hump and the densities 
given by (5.1). Note also the relationship 

N-l 

4 N h ~ 1) = Pn{x n+ i - X n ) . 
n=p 

As the number of cars increases (continuum limit) this leads to the formula 

lim e(N h - 1) = / p(x)dx = 7rp N + (1 - p N )A r , (5.2) 
where p = p(x) is as in (5.1) above and A r is the area under the function r = r(x). 

We now describe the behavior in the continuum limit of the problem solved by this script, using 
the theory of shocks and characteristics developed in section 2 here and elsewhere. The results of 
this analysis are built into the scheme graphics, that compare the actual solution of the equations 
(1.4) with the predictions here. The good agreement found is a confirmation of the correctness of 
the theory in section 2. In the continuum limit we use equation (2.5) to deal with the well behaved 
parts of the solution — where we can use the characteristic method — and equations (2.7) and 
(2.8) to deal with the discontinuities (shocks). 

Notice that in this case the wave speed satisfies c = 4 — Sp | and is a linear function of the density 



p. It then follows that in this example c is also a conserved quantity. Thus we can consider 
the solution of the continuum limit problem fully in terms of c. It is easy to see that c satisfies the 
equation 

dc d 1 

= ^ + ^(-c 2 ), with c{x,0) = c N -C{x), (5.3) 

where —4 < cn = 4 — 8pN < 4 and C = (4 + c/v)r(:r). Thus the initial profile for c has a "dip" 
instead of a "hump". In terms of c the shock condition (2.7) states: the shock speed is the 
average of the characteristic speeds on the sides of the shock. 



WARNING: this is true only for this case of a quadratic flow function Q = Q(p). 



Similarly, (2.8) becomes: across shocks c decreases — which is true for all flow functions Q. 



The characteristic curves are given by 



dx _ 
dt G 



dS _ _ c2 
dt ° 



on 



- with c constant. Furthermore 
them, where S = c x is the slope of the solution. This shows that S will eventually go to — oo on 
any characteristic where S starts negative, at time t = —1/S((, 0). Here £ is the value of x on the 
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dC 

characteristic at time t = and S(x,0) = — r~( x )- This follows from the general solution of the 

dx 

equation above for S along characteristics: 



S 



S(C,Q) 

l + 5(C,0)i 



An analysis of this problem shows that a shock will form — starting on the characteristic where 

dC 

S is negative and has the largest absolute value. Let this be given by S m = S(( m ,0) = — —(C)- 



dx 



Then the shock starts at 



ts = -q- and x s = (m + (c N - C(( m ))t s . 
Note that ( m must correspond to a location on the back end of the initial hump. 



Typical solution after shock forms. 

P 

Ps 




S N 
Areas under shock and multiple valued curve are equal. 



Typical solution after shock forms. 




* * 

*,* Areas enclosed by the curves are equal. 



Figure 5.2: Solution after a shock forms. The figure shows how the multiple values in the solution 
by characteristics are eliminated by the shock. The shock is located so that the area under the 
density curve (number of cars) is preserved. In this case, because the flow function is quadratic, 
the area under the wave velocity curve is also preserved. 



Eventually (time "large") most of the characteristics that start with x somewhere in the the position 
of the initial hump "die" at the shock. More precisely, the density disturbance will be made up 
only from characteristics that start near the leading edge (x = 0) of the initial hump, i.e. 

x = ct + C ) with c = Cat — C (() and £ small and negative. (5.4) 
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The mechanism behind this is simple. An inspection of the solution by characteristics (see figure 
5.2) shows that (as time advances) the initial hump in p (equivalently, in c) deforms. The back part 
steepens while the front part stretches. Eventually a shock forms on the back and all the details of 
the variations in the density are absorbed by it: the characteristics reach the shock and terminate 
there. The only part that remains is the very stretched out front. Because the stretching is linear 
in p, this part becomes a straight line, joining the front edge of the shock with the position of the 
leading characteristic starting at the front edge of the initial hump (i.e. x = c N t). Thus the wave 
takes a triangular form (backward saw-tooth) with the shock on the back and a corner moving 
at the characteristic speed cv at the front. Furthermore, because of conservation, the total 

area in the saw-tooth must be equal to that in the original hump. We can be a bit more precise 

x 

using (5.4), that yields the approximation c ~ — because £ is small. Then 

1. There is a shock at xs = c N t — \pttA. 

x 



(5.5) 



2. For xs < x < c N t c 

3. Elsewhere c = c N . 

4 - c 

4. The car density follows from p = — - — . 

8 

Here A is the area of the bump in c. That is 

A = / C(x)dx = (4 + c N )A r « 8(e(N h - 1) - irp N ) , 

where we have used (5.2) and the fact that 4 + c N = 8(1 — p N ) to write the last (approximate) 
equality. The formula for xs follows because this area must be conserved. Specifically, note that 
the value of c immediately ahead of the shock is given by 

x s 

cs = — = c N - 

The value of xs is selected so that the area of the triangular saw-tooth equals A. That is, so that 
the equation 2 A = (c N — cs)(c N t — xs) holds. Note that here we have used the fact that c 
itself is conserved, as stated earlier in (5.3). 

A more detailed justification of the arguments above can be found in the book by G. B. Whitham: 
Linear and Nonlinear Waves. Hopefully, it will also be done in the lectures. 




Simple Traffic Flow Model. MIT, Friday March 26, 1999 — Rosales. 19 

Of course, the triangular shape is achieved exactly only for t — >■ oo, so that these are not very good 
approximations as far as the script is concerned (which cannot run the solution for a very long 
time, particularly in the continuum limit when N has to be large). The main source of error occurs 
in terms of space and time translations. That is, the shape of the solution described by (5.5) is 
achieved fairly quickly, but it is not properly centered 10 . The source of these errors is that £ in (5.4) 
is small, but does not vanish. 

Thus, in order to do the graphical comparisons of the continuum limit with the actual solutions 
of the (1.4), the script quadCFSM uses an improved approximation, which we describe next 
(the idea is actually very simple). As stated earlier, after a while the details of the solution are 







\ x = x L ^^^^ jff I 


\ 

C = ( 


~ % \ / / 

\ / c ~ c.. + Bx 

Area = A \ / 


c 
- -4 



c = c(x, 0) and approximation near x = 0. 

Figure 5.3: The figure shows how the initial wave velocity profile 
is approximated, near the origin, by a straight line of slope B. 



determined only by the part of the initial data close to the origin, with the rest meaningful only 

as far as determining the area parameter A. Thus, instead of reducing the information about this 

10 Given any solution c = c(x, t), c = c(x — x*, t — i*) is also a solution (for any constants a;* and i*). The problem 
with the approximation (5.5) is mainly the existence of displacements a;* and it does not account for. 
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zone to just the fact that it is "near" x = where (at t = 0) c = cn, we approximate the initial 
data for x negative near the origin as follows (see figure 5.3): 

c(x, 0) ~ c N + Bx for some constant B > . (5.6) 

In terms of (5.3) this just means that we write C ~ — Bx. With initial conditions of this form, 



x 



equation (5.3) can be solved exactly. Thus, we replace the approximation c ~ — used earlier for the 
front part of the sawtooth, by the solution that follows from initial conditions as in (5.6). Other 
than this, we use the same ideas that lead to (5.5), to obtain the improved approximation: 



1. There is a shock at xs = c N t 

B 



2. For xs < x < c N t c = c N + 

where x x = c N t x and = —1/B. 



1 + Bt 



(x — c N t) 



I2A(1 + Bt) 



3. Elsewhere c = c N . 

4 - c 

4. The car density follows from p = — - — . 

8 



/ 2AB 1 / 2AB 

5. Right ahead of the shock: eg = c N — \ — and ps = Pn + „ 1 



1 + Bt 



8 V 1 + Bt 



(5.7) 



Just one issue remains now and it is how to best choose B. For a given target time around which one 
desires the approximation to be good, one can use (5.5) to get an estimate of what is the range of 
characteristics that are making up the front of the saw-tooth. That is, one can determine an interval 
xi < x < such that the characteristics originating there at time t = do not die at the shock up 
until after the target time. It is then over this range Xi < x < that we need the approximation 
(5.6) to hold. Of course, the target time cannot be too small, for the range xi < x < has to be 
small enough that an approximation of the form (5.6) makes sense. Once xi is determined, we can 
choose B by conservation of area (i.e. cars). That is, we require that the area under the straight 

line (5.6) be the same as the area under the curve it replaces. This yields the equation 

1 r° r° 

-Bx\= \ (cn — c(x,0))dx = 8 / (p(x,0) — pn)cIx . 

2, Jxi Jxi 

Now we can associate xi with the initial position of one of the cars in (1.4), say car number £, and 
make the following approximation 

, N-l 

(p(x, 0) - p N )dx w Pn{x n +\ ~ x n ) + x t p N = e(N - £) + x t p N . 



f 
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Thus, we obtain 

B = l6 e{N - £ \ +XipN . (5.8) 



Exact and approximate solutions. 





»\ 

VV 
<\ 


n 

r 


Shock 

\ 










Dashed line = approximation. 



X 



Figure 5.4: This figure shows an example comparing an exact so- 
lution of the continuum limit and the approximation given by (5.7). 



Remark 5.2 As final point: quadCFSM puts cars only a finite distance behind the initial hump; 
i.e. the initial conditions in (5.1) are defined only for rci(O) < x < x N (0) = 0, where Xi(0) < —ir is 
actually not too large in size 11 . Thus, if the computation were to be run for long enough, all the 
cars would eventually go through the disturbance and emerge in front of it where p is identically p^. 
That is, the solution of the ODE settles down to p n = p N after a sufficiently long time. This follows 
from the fact that the car speed u is always bigger than the wave speed c. 



11 The size of xi(0) is not fixed and is calculated to have enough cars in the problem to allow the saw-tooth shape 
enough time to develop fully before all the cars go through it. 



